This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
library(spatstat)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
install.packages("tidyverse")
library(tmap)
library(sf)
library(geojsonio)
library(tmaptools)
library(tidyverse)
library(rgdal)
library(raster)
library(SciViews)
library(base)
library(broom)
library(spdep)
####Basemap
library(tmap)
BoroughMap <- geojson_read("https://raw.githubusercontent.com/ft-interactive/geo-data/master/uk/london/london-wards-2014.geojson", what = "sp")
BNG = "+init=epsg:27700"
BoroughMapBNG <- spTransform(BoroughMap,BNG)
WGS = "+init=epsg:4326"
tmap_mode("view")
tm_shape(BoroughMapBNG) +tm_polygons(col = NA, alpha = 0.5)
#aggregate city of london
library(sf)
library(dplyr)
library(rgdal)
BoroughMapSF <- st_as_sf(BoroughMapBNG)
library(rgeos)
CityofLondon<- gUnionCascaded(BoroughMapBNG[c(630:654), ])
plot(CityofLondon)
CityofLondonSF <- st_as_sf(CityofLondon)
CityofLondonSF$gss_code_ward<-"E05000001" ##add same columns in city of london, for aggregating with other borough
CityofLondonSF$gss_code_borough <- "E09000001"
CityofLondonSF$borough <-"City Of London"
CityofLondonSF$ward <- "City Of London"
BoroughMapSFnew <- BoroughMapSF[-c(630:654), ]
LondonWardSF <- rbind(BoroughMapSFnew,CityofLondonSF) ##Attention New not new
LondonWard <- as(LondonWardSF,"Spatial")
rownames(LondonWardSF) <- 1:nrow(LondonWardSF)
plot(LondonWard)
tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)
####population density (using https://data.london.gov.uk/dataset/land-area-and-population-density-ward-and-borough 2016)
library(tidyverse)
Density <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/population_density_2016_final.csv",na = "n/a")
#normalization population_per_square_km
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
density.normalization <- normalization(Density$`population density`)
Density$Density.normalization <- density.normalization
#join the data to basemap P3P27
LondonWard.join.popdens <- LondonWardSF%>% left_join(Density,by=c("gss_code_ward"="New Code"))
qtm(LondonWard.join.popdens,fill="Density.normalization")
#### coexistence of old and new building
#read building list
building <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Listed Buildings/ListedBuildings_06Dec2018.shp",layer="ListedBuildings_06Dec2018")
library(maptools)
LondonBoundary<-unionSpatialPolygons(LondonWard,gss_code_borough)
plot(LondonBoundary)
### Public transport accessibility level(PTAL)
PTAL <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Ward2014 Avpublic transport accessibility level2015.csv",na = "n/a")
#normalization PTAL
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
PTAL.normalization <- normalization(PTAL$AvPTAI2015)
PTAL$PTAL.normalization <- PTAL.normalization
#join the data to basemap P3P27
LondonWard.join.PTAL <- LondonWardSF%>% left_join(PTAL,by=c("gss_code_ward"="Ward Code"))
qtm(LondonWard.join.PTAL,fill="PTAL.normalization")
qtm(LondonWard.join.PTAL,fill="AvPTAI2015")
###Public open space accessibility level(POSAL)
POSAL <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/access_public_open_space_ward.csv",na = "n/a")
POSAL$POSAL.normalization <- normalization(POSAL$`Open Space`)
LondonWard.join.POSAL <- LondonWardSF%>% left_join(POSAL,by=c("gss_code_ward"="WD13CD"))
#qtm(LondonWard.join.POSAL,fill="POSAL.normalization")
#combine two accessibility sub-index
Access.normalization<- merge(POSAL,PTAL, by.x="WD13CD",by.y="Ward Code",all=TRUE)
#remove useless column and rows
Access.normalization[,c(4:8,10:11)] <- NULL
#calculate the mean of two accessibility level
Access.normalization$access.mean <- rowMeans(Access.normalization[c('PTAL.normalization', 'POSAL.normalization')], na.rm=TRUE)
economic vitality
###economic vitality
#CLASS 1 house price
#HousePrice <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/house_price_mean2016.csv",na = "n/a")
#HousePrice.normalization <- normalization(HousePrice$Value)
#HousePrice$HousePrice.normalization <- HousePrice.normalization
#CLASS 2 economic activity
EcoAct <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/economic_activity_ward.csv",na = "n/a")
EcoAct.normalization <- normalization(EcoAct$`Economically active: Total`)
EcoAct$EcoAct.normalization <- EcoAct.normalization
land use mix
### Land use mix
#https://en.wikipedia.org/wiki/Planning_use_classes_in_England
# six classes--
landcover <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/gis_osm_landuse_a_free_1.shp",layer="gis_osm_landuse_a_free_1")
landcoverBNG <- spTransform(landcover, BNG)
landcoverBNGSF <- st_as_sf(landcoverBNG)
AOIS <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/gis_osm_pois_a_free_1.shp",layer="gis_osm_pois_a_free_1")
AOISBNG <- spTransform(AOIS, BNG)
AOISBNGSF <- st_as_sf(AOISBNG)
#find missing data in each column in landcoverBNGSF
colSums(is.na(landcoverBNGSF))
#get the classes of landcover
summary(landcoverBNGSF)
#CLASS1 greenspace
greenspaceBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="park" | landcoverBNGSF$fclass =="forest"|landcoverBNGSF$fclass =="nature_reserve"|landcoverBNGSF$fclass =="grass"),]
greenspaceBNGSP <- as(greenspaceBNGSF,"Spatial")
#CLASS2 industrial and commercial
industrialBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="industrial"| landcoverBNGSF$fclass =="commercial"),]
industrialBNGSP <- as(industrialBNGSF,"Spatial")
#CLASS 3 residential
residentialBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="residential" ),]
residentialBNGSP <- as(residentialBNGSF,"Spatial")
#CLASS 4 recreation-leisure
recreationBNGSF1 <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="recreation_ground"| landcoverBNGSF$fclass =="allotments"),]
recreationBNGSF2 <- AOISBNGSF[ which(AOISBNGSF$fclass =="attraction"|AOISBNGSF$fclass =="theatre"|AOISBNGSF$fclass =="arts_centre"),]
recreationBNGSF <- rbind(recreationBNGSF2,recreationBNGSF1)
recreationBNGSP <- as(recreationBNGSF,"Spatial")
#CLASS 5 retail
retailBNGSF1 <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="retail" ),]
retailBNGSF2 <- AOISBNGSF[ which( AOISBNGSF$fclass =="retail" ),]
retialBNGSF <- rbind(retailBNGSF2,retailBNGSF1)
retailBNGSP <- as(retailBNGSF,"Spatial")
#CLASS 6--community services
communityBNGSF <- AOISBNGSF[ which(AOISBNGSF$fclass =="school"|AOISBNGSF$fclass =="hospital"|AOISBNGSF$fclass =="nursing_home"|AOISBNGSF$fclass =="university" ),]
communityBNGSP <- as(communityBNGSF,"Spatial")
#CLASS 7--transport
transport <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/road_polygon_merge.shp",layer="road_polygon_merge")
transportBNG <- spTransform(transport, BNG)
transportBNG <- transportBNG[transportBNG@data$AREA >10000,]
# Assign the land use area for each subset. Here we first use the most widespread residential as example
library(raster)
r <- raster(ncols=400, nrows=400) #generate raster size(quantity)
extent(r) <- extent(residentialBNGSP) #Important! Assign the extent of raster to cover the same extents of the polygon
residential.Raster <- rasterize(residentialBNGSP,r,background=NA) #convert polygon to raster
#tmap_mode("view")
#qtm(residential.Raster)+tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)
#extract rasters from polygons#(extract() cannot use method=bilinear in raster-in-polygon)
residential.extract <- raster::extract(residential.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
#remove rows with NA
row.has.na <- apply(residential.extract, 1, function(x){any(is.na(x))})
residential.extract2 <-residential.extract[!row.has.na,]
#merge the pixel in same ward(aggregate the frequency of the same data in one)
library(plyr)
residential.extract.sum <- factor(residential.extract2$ID)
residential.extract.sum <- table(residential.extract2$ID)
residential.extract.wardsum <- as.data.frame(residential.extract.sum)
#create new order number for LondonWardSF, for merge the residential.extract and LondonWardSF by order
LondonWard2SF <- LondonWardSF
LondonWard2SF$Var1 <- 1:630
LondonWard2SF[,'Var1']<-factor(LondonWard2SF[,'Var1'])
#merge two dataframe
residential.extract.final <- merge(residential.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
#residential.extract.final [!duplicated(residential.extract.final ), ]
#replace NA values with 0
residential.extract.final[is.na(residential.extract.final)] <- 0
#rename the column
colnames(residential.extract.final)[2] <- "Freq.residential"
#raster cell size
#we have to convert raster to polygon because of raster(). This function only compute the pixel area in longitude/latitude coordiante system
residential.repolygon <- rasterToPolygons(residential.Raster, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=TRUE)
#qtm(residential.repolygon)
residential.repolygonSF <- st_as_sf(residential.repolygon)
residential.repolygonSF$area <- st_area(residential.repolygonSF)
#so the single pixel area is 28999.54 [m^2](0.029 km^2)
land use mix functions
##based on the details above, we create functions for land use analysis
function.raster <- function(landsp){
extent(r) <- extent(landsp)
landsp.Raster <- rasterize(landsp,r,background=NA)
return(landsp.Raster)
}
function.extract<- function(landsp.Raster){
landsp.extract <- raster::extract(landsp.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
row.has.na <- apply(landsp.extract, 1, function(x){any(is.na(x))})
landsp.extract2 <-landsp.extract[!row.has.na,]
landsp.extract.sum <- factor(landsp.extract2$ID)
landsp.extract.sum <- table(landsp.extract2$ID)
landsp.extract.wardsum <- as.data.frame(landsp.extract.sum)
landsp.extract.final <- merge(landsp.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
landsp.extract.final[is.na(landsp.extract.final)] <- 0
return(landsp.extract.final)
}
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(index.dataframeSP )+
tm_polygons(c("Index", "MixIndex","Density.normalization","EcoAct.normalization","access.mean")) +
tm_facets(sync = TRUE, ncol = 1,nrow=5,drop.empty.facets = F)
The number of facets exceeds the limit of 4. The limit can be extended to 5 with:
tmap_options(limits = c(facets.view = 5))